
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Empirical Application(Timber Dataset)               %
% Test: Constancy of the Slope Coefficients           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all

load('data.mat'); % load timber dataset
vol2=vol2./1000; 
vol3=vol3./1000;
apv2=apv2./100;
apv3=apv3./100;

L2=size(wb2,1); % number of auctions with N=2
L3=size(wb3,1); % number of auctions with N=3
L = L2+L3;

% covariates
x2=[ones(L2,1),apv2,vol2]; 
x3=[ones(L3,1),apv3,vol3];
x = [x2;x3];
p=size(x,2);

% winning bids
wb=[wb2;wb3];

N=[2,3]; % Number of bidders
alpha=0.12:0.02:0.80; % quantile levels alpha
psi2=ones(L2,1)*(N(1)*(alpha.^(N(1)-1))-(N(1)-1)*(alpha.^N(1))); % Psi function given N=2
psi3=ones(L3,1)*(N(2)*(alpha.^(N(2)-1))-(N(2)-1)*(alpha.^N(2))); % Psi function given N=3
psi = [psi2;psi3];
m=size(alpha,2);

B=5000; % Bootstrap replications
toler=2; % Parameter of tolerance for precision of the estimation

gamma0=ones(1,m); % Initial guess for the intercept
gammap=ones(p-1,1); % Initial guess for the slope
gamma=ones(p,1);

% Computation of the Constancy of the Slope Coefficients Test Statistic

% Estimation under H0
[obj_cqr,g0,gp,res_cqr]=CQR(psi, p, wb, x, toler,gamma0,gammap); 
obj0=sum(obj_cqr)/m; % objective function under H0
gamma0=g0; % estimate of the intercept under H0
gammap=gp; % estimate of the slope under H0

% Estimation under H1
for k=1:m
    [obj,res,~]=MM(psi(:,k), p, wb, x, toler, gamma);
    resk(:,k)=res; % residuals under H1
    objk(k,:)=obj; % objective function for each k under H1
end
    obj1=sum(objk)/m;
    M=obj0-obj1; % CQR test statistic

% Random Weighting Bootstrap of the Test Statistic for Constancy of the Slope Coefficients 

for b=1:B
        
        w = [mnrnd(L2,ones(L2,1)/L2)';mnrnd(L3,ones(L3,1)/L3)']; % weights
         
        [obj_cqr_wb,~,~,~]=CQR_w(psi, p, wb, x, w, toler,gamma0,gammap);
        obj0wb=sum(obj_cqr_wb)/m; % random weighted objective function under H0
        
   for k=1:m  
         
        res_cqrw=res_cqr(:,k).*w; % residuals under H0 weighted by w
        resw=resk(:,k).*w; % residuals under H1 weighted by w    
        obj_cqrw(k,:)=sum(psi(:,k).*res_cqrw)-sum(res_cqrw(res_cqrw<0)); % objective function under H0 weighted by w
        objw(k,:)=sum(psi(:,k).*resw)-sum(resw(resw<0)); % objective function under H1 weighted by w
        
        % objective functions under the random weighting method
        [obj_wb,~,~]=MMw(psi(:,k), p, wb, x, w, toler, gamma);
        objk_wb(k,:)=obj_wb;
   end
        obj0w=sum(obj_cqrw)/m;
        obj1w=sum(objw)/m;
        Mw=obj0w-obj1w;
        
        obj1_wb=sum(objk_wb)/m;
        Ms=obj0wb-obj1_wb;
        Mstar(:,b)=Ms-Mw; % Bootstraped test statistic
        b
end

% Critical Values
c99=quantile(Mstar',0.99);
c95=quantile(Mstar',0.95);
c90=quantile(Mstar',0.90);
    
% P-value computation    
R=zeros(1,B);
for b=1:B
    if Mstar(:,b)>=M
            R(:,b)=1;
    else
            R(:,b)=0;
    end
end
pvalue=sum(R')/B
    


    
    